% csv_twoD_extractor.m  3/20/2014  Parker MacCready
%
% this is for extracting the results of ROMS simulations
%
% it is focused on time series of 2D fields
%
% it saves things in a csv file

clear;addpath('../alpha'); [Tdir] = toolstart;

% &&&&&&&&&&& USER EDIT THIS &&&&&&&&&&&&&&
doLocal = true;
if doLocal
    indir = '/Users/PM3/Documents/roms/output/';
else
    indir = '/pmraid1/daves/runs/'; % for skua
end
basename = 'salish_2006_4';
% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

%% create top level output directory if needed
outdir = [Tdir.output,'twoD_out/'];
if ~exist(outdir,'dir'); mkdir(outdir); end;

%% choose run to work on and set the basename
disp('********* twoD_extractor.m ****************')
disp(' ')
if doLocal
    if 0 % interactive version
        disp('** Select folder where history files are **');
        pth = uigetdir(indir);
    else % hardwired
        pth = [indir,basename,'/'];
    end
else % hardwired for skua
    pth = [indir,basename,'/OUT/'];
end

% get the run definition
run = roms_createRunDef('my_run',pth);

% &&&&&& USER EDIT IF WANTED &&&&&&
if 1 % a limited time span
    year = datestr(run.his.nctime(1),'yyyy');
    year = str2num(year);
    t0 = datenum(year,6,18);
    t1 = datenum(year,6,19);
else % the full length of the run
    t0 = run.his.nctime(1);
    t1 = run.his.nctime(end);
end

%% prepare output arrays

lon = run.grid.lon;
lat = run.grid.lat;
mask = run.grid.mask;
lonu = run.grid.lonu;
latu = run.grid.latu;
lonv = run.grid.lonv;
latv = run.grid.latv;
[nr,nc] = size(lon);

flist = {'salt_bot','temp_bot','salt_top','temp_top','u_bot','v_bot', ...
    'u_bar','v_bar','u_top','v_top'};
plist = {'max','min','mean','std'};
zmat = zeros(nr,nc);
nmat = nan(nr,nc);

% initialize all as zeros or nans
for ff = 1:length(flist)
    for pp = 1:length(plist)
        ppp = plist{pp};
        fff = flist{ff};
        fname = [fff,'_',ppp];
        switch ppp
            case {'max','min'}
                eval([fname,' = nmat;']);
            case {'mean','std'}
                eval([fname,' = zmat;']);
        end
    end
end

%% do the extraction and averaging

it0 = dsearchn(run.his.nctime',t0);
it1 = dsearchn(run.his.nctime',t1);
ncn = run.his.ncn(it0:it1);
NT = length(ncn);
NZ = length(run.grid.cs);

%% first find the means

plist1 = {'mean'};
disp(' --- working on means ---')
for tt = ncn(1):ncn(end)
    
    % name the fn
    tt_str = ['0000',num2str(tt)];
    tt_str = tt_str(end-3:end);
    fn = [run.his.dirname,run.his.basename,tt_str,'.nc'];
    
    % some things we only need once
    if tt == ncn(1)
        H = run.grid.H;
        S = Z_get_S(fn);
        z = Z_s2z_rho(H,0*H,S);
        mab = squeeze(z(1,:,:)) + H; % mab of the deepest rho point
    end
    
    % screen output about progress
    if mod(tt,5)==0; disp(['  fn = ',fn]); end;
    
    % get the raw fields
    salt_bot = nc_varget(fn,'salt',[0 0 0 0],[1 1 -1 -1]);
    temp_bot = nc_varget(fn,'temp',[0 0 0 0],[1 1 -1 -1]);
    salt_top = nc_varget(fn,'salt',[0 NZ-1 0 0],[1 1 -1 -1]);
    temp_top = nc_varget(fn,'temp',[0 NZ-1 0 0],[1 1 -1 -1]);
    ur_bot = nc_varget(fn,'u',[0 0 0 0],[1 1 -1 -1]);
    vr_bot = nc_varget(fn,'v',[0 0 0 0],[1 1 -1 -1]);
    ur_bar = nc_varget(fn,'ubar');
    vr_bar = nc_varget(fn,'vbar');
    ur_top = nc_varget(fn,'u',[0 NZ-1 0 0],[1 1 -1 -1]);
    vr_top = nc_varget(fn,'v',[0 NZ-1 0 0],[1 1 -1 -1]);
    
    % set everything on land to zero
    % (only important for velocity, when we interpolate)
    rlist = {'ur_bot','vr_bot','ur_bar','vr_bar','ur_top','vr_top'};
    for ff = 1:length(rlist)
        fff = rlist{ff};
        eval([fff,'(isnan(',fff,')) = 0;']);
    end
    
    % interpolate to the rho grid
    u_bot = interp2(lonu,latu,ur_bot,lon,lat);
    v_bot = interp2(lonv,latv,vr_bot,lon,lat);
    u_bar = interp2(lonu,latu,ur_bar,lon,lat);
    v_bar = interp2(lonv,latv,vr_bar,lon,lat);
    u_top = interp2(lonu,latu,ur_top,lon,lat);
    v_top = interp2(lonv,latv,vr_top,lon,lat);
    
    for ff = 1:length(flist)
        for pp = 1:length(plist1)
            ppp = plist1{pp};
            fff = flist{ff};
            fname = [fff,'_',ppp];
            switch ppp
                case 'mean'
                    eval([fname,'=',fname,'+',fff,'/NT;']);
            end
        end
    end
end

%% now find all the rest of the items on plist

plist2 = {'max','min','std'};
disp(' --- working on max, min, std ---')
for tt = ncn(1):ncn(end)
    
    % name the fn
    tt_str = ['0000',num2str(tt)];
    tt_str = tt_str(end-3:end);
    fn = [run.his.dirname,run.his.basename,tt_str,'.nc'];
    
    % screen output about progress
    if mod(tt,5)==0; disp(['  fn = ',fn]); end;
    
    % get the raw fields
    salt_bot = nc_varget(fn,'salt',[0 0 0 0],[1 1 -1 -1]);
    temp_bot = nc_varget(fn,'temp',[0 0 0 0],[1 1 -1 -1]);
    salt_top = nc_varget(fn,'salt',[0 NZ-1 0 0],[1 1 -1 -1]);
    temp_top = nc_varget(fn,'temp',[0 NZ-1 0 0],[1 1 -1 -1]);
    ur_bot = nc_varget(fn,'u',[0 0 0 0],[1 1 -1 -1]);
    vr_bot = nc_varget(fn,'v',[0 0 0 0],[1 1 -1 -1]);
    ur_bar = nc_varget(fn,'ubar');
    vr_bar = nc_varget(fn,'vbar');
    ur_top = nc_varget(fn,'u',[0 NZ-1 0 0],[1 1 -1 -1]);
    vr_top = nc_varget(fn,'v',[0 NZ-1 0 0],[1 1 -1 -1]);
    
    % set everything on land to zero
    % (only important for velocity, when we interpolate)
    rlist = {'ur_bot','vr_bot','ur_bar','vr_bar','ur_top','vr_top'};
    for ff = 1:length(rlist)
        fff = rlist{ff};
        eval([fff,'(isnan(',fff,')) = 0;']);
    end
    
    % interpolate to the rho grid
    u_bot = interp2(lonu,latu,ur_bot,lon,lat);
    v_bot = interp2(lonv,latv,vr_bot,lon,lat);
    u_bar = interp2(lonu,latu,ur_bar,lon,lat);
    v_bar = interp2(lonv,latv,vr_bar,lon,lat);
    u_top = interp2(lonu,latu,ur_top,lon,lat);
    v_top = interp2(lonv,latv,vr_top,lon,lat);
    
    for ff = 1:length(flist)
        for pp = 1:length(plist2)
            ppp = plist2{pp};
            fff = flist{ff};
            fname = [fff,'_',ppp];
            switch ppp
                case 'max'
                    eval([fname,'=max(',fname,',',fff,');']);
                case 'min'
                    eval([fname,'=min(',fname,',',fff,');']);
                case 'std'
                    eval([fname,'=',fname,'+(',fff,'-',fff,'_mean).^2/(NT-1);']);
            end
        end
    end
end

%% take square roots of std
plist3 = {'std'};
for ff = 1:length(flist)
    for pp = 1:length(plist3)
        ppp = plist3{pp};
        fff = flist{ff};
        fname = [fff,'_',ppp];
        eval([fname,' = sqrt(',fname,');']);
    end
end

%% mask out land
for ff = 1:length(flist)
    for pp = 1:length(plist)
        ppp = plist{pp};
        fff = flist{ff};
        fname = [fff,'_',ppp];
        eval([fname,'(~mask) = NaN;']);
    end
end

%% test plotting
if 0
    close all; Z_fig(10);
    
    if 0
        flist9 = flist;
        plist9 = plist;
    else
        flist9 = {'salt_top','temp_top','u_top','v_top'};
        plist9 = {'max','min','mean','std'};
    end
    
    count = 1;
    for ff = 1:length(flist9)
        for pp = 1:length(plist9)
            ppp = plist9{pp};
            fff = flist9{ff};
            fname = [fff,'_',ppp];
            eval(['varb = ',fname,';']);
            disp([' === plotting ',fname]);
            subplot(length(plist9),length(flist9),count)
            Z_pcolorcen(lon,lat,varb);
            hold on
            Z_dar
            Z_addcoast('combined',Tdir.coast);
            title(strrep(fname,'_',' '));
            colorbar
            count = count + 1;
        end
    end
end

%% save output files

% create output directory if needed
outdir2 = [outdir,basename,'/'];
if ~exist(outdir2,'dir'); mkdir(outdir2); end;

if 0
    flist10 = flist;
    plist10 = plist;
else
    flist10 = {'salt_top'};
    plist10 = {'max'};
end

for ff = 1:length(flist10)
    for pp = 1:length(plist10)
        ppp = plist10{pp};
        fff = flist10{ff};
        fname = [fff,'_',ppp];
        eval(['varb = ',fname,';']);
        outname = [outdir,basename,'/',fname,'.csv'];
        fid = fopen(outname,'w+');
        for rr = 1:nr
            for cc = 1:nc
                if ~isnan(varb(rr,cc))
                    xx = lon(rr,cc);
                    yy = lat(rr,cc);
                    val = varb(rr,cc);
                    fprintf(fid,'%0.5f,%0.5f,%0.5f\n',xx,yy,val);
                end
            end
        end
        fclose(fid);
    end
end

%% and finally write out the bathymetry and mab in the same way
fname_list = {'H','mab'};
for ff = 1:length(fname_list)
    fname = fname_list{ff};
    eval(['varb = ',fname,';']);
    outname = [outdir,basename,'/',fname,'.csv'];
    fid = fopen(outname,'w+');
    for rr = 1:nr
        for cc = 1:nc
            if mask(rr,cc)
                xx = lon(rr,cc);
                yy = lat(rr,cc);
                val = varb(rr,cc);
                fprintf(fid,'%0.5f,%0.5f,%0.5f\n',xx,yy,val);
            end
        end
    end
    fclose(fid);
end



